---
title: "PSRM Replication for Analyze the Attentive and Bypass Bias"
output: html_document
---

# Load packages

```{r}
library(haven)
library(estimatr)
library(margins)
library(tidyverse)
library(stargazer)
library(rmeta)
library(hrbrthemes)
library(modelsummary)
```

# Load data

```{r}
load("R_datasets.RData")
```

# Table 4 -----------------------------------------------------------------

```{r}
m0 <- lm(outcome ~ treatment*mvscore, lucid_long)
m0.robust <- commarobust(m0, se_type = "CR2", clusters = 
                          na.omit(lucid_long[,c("treatment", "outcome", "mvscore","subject")])$subject)

m0.robust %>%
  modelsummary()
```

# Figure 3 ----------------------------------------------------------------

```{r}
# get proportion of correct responses within each MV category
lucid_long %>% 
 group_by(MV, mvscore) %>% 
 summarize(size = n()) %>% 
 group_by(MV, add = TRUE) %>% 
 mutate(proportion = size/sum(size)) -> proportion

# calculate ATEs
map_df(list(lm(outcome ~ treatment + mvscore*treatment, lucid_long, MV == 1 & EXP == 1),
            lm(outcome ~ treatment + mvscore*treatment, lucid_long, MV == 1 & EXP == 2),
            lm(outcome ~ treatment + mvscore*treatment, lucid_long, MV == 1 & EXP == 3),
            lm(outcome ~ treatment + mvscore*treatment, lucid_long, MV == 1 & EXP == 4),
            lm(outcome ~ treatment + mvscore*treatment, lucid_long, MV == 2 & EXP == 1),
            lm(outcome ~ treatment + mvscore*treatment, lucid_long, MV == 2 & EXP == 2),
            lm(outcome ~ treatment + mvscore*treatment, lucid_long, MV == 2 & EXP == 3),
            lm(outcome ~ treatment + mvscore*treatment, lucid_long, MV == 2 & EXP == 4),
            lm(outcome ~ treatment + mvscore*treatment, lucid_long, MV == 3 & EXP == 1),
            lm(outcome ~ treatment + mvscore*treatment, lucid_long, MV == 3 & EXP == 2),
            lm(outcome ~ treatment + mvscore*treatment, lucid_long, MV == 3 & EXP == 3),
            lm(outcome ~ treatment + mvscore*treatment, lucid_long, MV == 3 & EXP == 4),
            lm(outcome ~ treatment + mvscore*treatment, lucid_long, MV == 4 & EXP == 1),
            lm(outcome ~ treatment + mvscore*treatment, lucid_long, MV == 4 & EXP == 2),
            lm(outcome ~ treatment + mvscore*treatment, lucid_long, MV == 4 & EXP == 3),
            lm(outcome ~ treatment + mvscore*treatment, lucid_long, MV == 4 & EXP == 4)), function(x)
             summary(margins(x, at = list(mvscore = c(0, 1, 2, 3))))[5:8,c("AME", "SE")], 
       .id = "model") %>% mutate(CORRECT = rep(0:3, 16), 
                                 MV = rep(1:4, each = 16), 
                                 EXP = rep(1:4, each = 4, length.out = 64)) -> results

# merge MV proportions and AMEs
merged <- merge(proportion, results, by.x = c("MV", "mvscore"),
                by.y = c("MV", "CORRECT"))

# labels
merged$EXP <- factor(merged$EXP, 
                     levels = 1:4, 
                     labels = c("Student Loans", 
                                "KKK", 
                                "Welfare", 
                                "Immigration"))

merged$MV <- factor(merged$MV, 
                    levels = 1:4, 
                    labels = c("MV: Scientific Publishing", 
                               "MV: Stadium Licenses", 
                               "MV: Sulfur Dioxide", 
                               "MV: Hazardous Plants"))

distance <- 0

# plots
merged %>% 
 ggplot(aes(x = mvscore, y = AME, shape = EXP, group = EXP)) + 
 geom_point(position = position_dodge(distance), size = 2, alpha = .85) + 
 geom_line(position = position_dodge(distance)) +
 geom_errorbar(aes(ymin = AME - 1.96*SE, ymax = AME + 1.96*SE), width = 0, 
               position = position_dodge(distance)) + 
 facet_wrap(. ~ MV, scale = "fixed") + 
 theme_minimal() + 
 geom_histogram(aes(x = mvscore, y = proportion), 
                stat = "identity", fill = "black", size = 0, alpha = .1) + 
 scale_y_continuous(sec.axis = sec_axis(~./4 * 100, 
                                        name = "Percent of Sample")) + 
 geom_hline(yintercept = 0) + 
 labs(x = "Number of Correct MVCs", 
      y = "Treatment Effect", 
      shape = "Experiment") + 
 theme(legend.position = "bottom") + guides(color = FALSE) + 
 scale_shape_manual(values = c(19, 17, 21, 24)) -> p

p + theme(axis.title.x = element_text(face = "bold"), 
          axis.title.y = element_text(face = "bold"),
          strip.text.x = element_text(face = "bold"))
```

# Figure C1 ---------------------------------------------------------------

```{r}
tibble(group = as.factor(c("MVC Non-Passers", "MVC Passers",
                           "MVC Non-Passers", "MVC Passers", 
                           "MTurk DiD Estimate", 
                           "Qualtrics DiD Estimate")),
       study = c("MTurk 1 Treatment Effects", 
                 "MTurk 1 Treatment Effects", 
                 "Qualtrics Treatment Effects", 
                 "Qualtrics Treatment Effects", 
                 "Difference-in-Differences Estimate", 
                 "Difference-in-Differences Estimate"),
       type = as.factor(c(rep("  ", 2), c(" ", " "), 
                          rep("", 1), "")),
       estimate = c(-.41, -.63, 
                    -.36, -.51, 
                    abs(c(-.31, -.15))),
       ses = c(.28, .17, .22, .15, .18, .14)) %>%
 mutate(group = forcats::fct_rev(group)) %>%
 ggplot(aes(x = group, y = estimate, 
            group = study)) +
 facet_wrap(. ~ type + study, nrow = 4, scales = "free_y") + 
 geom_pointrange(aes(ymin = estimate - 1.96*ses, ymax = estimate + 1.96*ses),
                 position = position_dodge(.5)) +
 coord_flip() +
 theme_minimal() +
 labs(y = "Estimate", x= "") +
 geom_hline(yintercept = 0) +
 theme(legend.position = "", strip.text = element_text(hjust = 0)) +
 ggtitle("Treatment Effects and DiD Estimates")
```

# Table E2 ----------------------------------------------------------------

```{r}
# code duration
lucid_wide$duration <- as.numeric(lucid_wide$V4 - lucid_wide$V3)

# code MVC scores for the 2nd round MV
mvscore <- NA
mvscore[which(lucid_wide$MVround2 == 1)] <- lucid_wide$mv1scale[which(lucid_wide$MVround2 == 1)]
mvscore[which(lucid_wide$MVround2 == 2)] <- lucid_wide$mv2scale[which(lucid_wide$MVround2 == 2)]
mvscore[which(lucid_wide$MVround2 == 3)] <- lucid_wide$mv3scale[which(lucid_wide$MVround2 == 3)]
mvscore[which(lucid_wide$MVround2 == 4)] <- lucid_wide$mv4scale[which(lucid_wide$MVround2 == 4)]

# min-max coding
mvscore <- (mvscore - min(mvscore, na.rm = T))/
 (max(mvscore, na.rm = T) - min(mvscore, na.rm = T))

# select outcome variables
e1_control <- lucid_wide$Q144_3; e1_treatment <- lucid_wide$Q146_3
e1_outcome <- lucid_wide$Q148_3; e1_survey <- lucid_wide$duration
e1_fmc <- lucid_wide$fmc1

e2_control <- lucid_wide$Q184_3; e2_treatment <- lucid_wide$Q182_3
e2_outcome <- lucid_wide$Q186_3; e2_survey <- lucid_wide$duration
e2_fmc <- lucid_wide$fmc2

e3_control <- lucid_wide$Q54_3; e3_treatment <- lucid_wide$Q12_3
e3_outcome <- lucid_wide$Q11_3; e3_survey <- lucid_wide$duration
e3_fmc <- lucid_wide$fmc3

e4_control <- lucid_wide$Q160_3; e4_treatment <- lucid_wide$Q162_3
e4_outcome <- lucid_wide$Q164_3 + lucid_wide$Q166_3 + lucid_wide$Q168_3
e4_survey <- lucid_wide$duration
e4_fmc <- lucid_wide$fmc4

# estimate models
e1_c <- lm(log(e1_control) ~ mvscore + factor(lucid_wide$MVround2), 
           subset = lucid_wide$EXPround2 == 1)
e1_t <- lm(log(e1_treatment) ~ mvscore + factor(lucid_wide$MVround2), 
           subset = lucid_wide$EXPround2 == 1)
e1_o <- lm(log(e1_outcome) ~ mvscore + factor(lucid_wide$MVround2), 
           subset = lucid_wide$EXPround2 == 1)
e1_s <- lm(log(e1_survey) ~ mvscore + factor(lucid_wide$MVround2), 
           subset = lucid_wide$EXPround2 == 1)
e1_f <- lm(e1_fmc ~ mvscore + factor(lucid_wide$MVround2), 
           subset = lucid_wide$EXPround2 == 1)
e1_m <- lm(log(lucid_wide$Q10_3 + .01) ~ mvscore + factor(lucid_wide$MVround2), 
           subset = lucid_wide$EXPround2 == 1)

e2_c <- lm(log(e2_control) ~ mvscore + factor(lucid_wide$MVround2), 
           subset = lucid_wide$EXPround2 == 2)
e2_t <- lm(log(e2_treatment) ~ mvscore + factor(lucid_wide$MVround2), 
           subset = lucid_wide$EXPround2 == 2)
e2_o <- lm(log(e2_outcome) ~ mvscore + factor(lucid_wide$MVround2), 
           subset = lucid_wide$EXPround2 == 2)
e2_s <- lm(log(e2_survey) ~ mvscore + factor(lucid_wide$MVround2), 
           subset = lucid_wide$EXPround2 == 2)
e2_f <- lm(e2_fmc ~ mvscore + factor(lucid_wide$MVround2), 
           subset = lucid_wide$EXPround2 == 2)
e2_m <- lm(log(lucid_wide$Q125_3) ~ mvscore + factor(lucid_wide$MVround2), 
           subset = lucid_wide$EXPround2 == 2)

e3_c <- lm(log(e3_control) ~ mvscore + factor(lucid_wide$MVround2), 
           subset = lucid_wide$EXPround2 == 3)
e3_t <- lm(log(e3_treatment) ~ mvscore + factor(lucid_wide$MVround2), 
           subset = lucid_wide$EXPround2 == 3)
e3_o <- lm(log(e3_outcome) ~ mvscore + factor(lucid_wide$MVround2), 
           subset = lucid_wide$EXPround2 == 3)
e3_s <- lm(log(e3_survey) ~ mvscore + factor(lucid_wide$MVround2), 
           subset = lucid_wide$EXPround2 == 3)
e3_f <- lm(e3_fmc ~ mvscore + factor(lucid_wide$MVround2), 
           subset = lucid_wide$EXPround2 == 3)
e3_m <- lm(log(lucid_wide$Q133_3) ~ mvscore + factor(lucid_wide$MVround2), 
           subset = lucid_wide$EXPround2 == 3)

e4_c <- lm(log(e4_control) ~ mvscore + factor(lucid_wide$MVround2), 
           subset = lucid_wide$EXPround2 == 4)
e4_t <- lm(log(e4_treatment) ~ mvscore + factor(lucid_wide$MVround2), 
           subset = lucid_wide$EXPround2 == 4)
e4_o <- lm(log(e4_outcome) ~ mvscore + factor(lucid_wide$MVround2), 
           subset = lucid_wide$EXPround2 == 4)
e4_s <-lm(log(e4_survey) ~ mvscore + factor(lucid_wide$MVround2), 
          subset = lucid_wide$EXPround2 == 4)
e4_f <-lm(e4_fmc ~ mvscore + factor(lucid_wide$MVround2), 
          subset = lucid_wide$EXPround2 == 4)
e4_m <-lm(log(lucid_wide$Q152_3) ~ mvscore + factor(lucid_wide$MVround2), 
          subset = lucid_wide$EXPround2 == 4)

time_results <- exp(bind_rows(c(coef(e1_m)[2], coef(e1_c)[2], 
                                coef(e1_t)[2], coef(e1_o)[2], coef(e1_s)[2]),
                              c(coef(e2_m)[2], coef(e2_c)[2], coef(e2_t)[2], 
                                coef(e2_o)[2], coef(e2_s)[2]),
                              c(coef(e3_m)[2], coef(e3_c)[2], coef(e3_t)[2], 
                                coef(e3_o)[2], coef(e3_s)[2]),
                              c(coef(e4_m)[2], coef(e4_c)[2], coef(e4_t)[2], 
                                coef(e4_o)[2], coef(e4_s)[2])) -1)*100 

pr_pass_fmc <- c(coef(e1_f)[2], coef(e2_f)[2], 
                 coef(e3_f)[2], coef(e4_f)[2])

time_results %>%
 bind_cols(pr_pass_fmc)
```

# Table F1 ----------------------------------------------------------------

```{r}
lucid_long %>%
  group_by(EXP) %>%
  do(mod = lm_robust(outcome ~ treatment + mvscore +
                       treatment*mvscore + factor(MV) + 
                       factor(round), ., clusters = subject)) -> model_list

modelsummary(model_list$mod)
```

# Table F2 -----------------------------------------------------------------

```{r}
m1 <- lm(outcome ~ treatment*mvscore*factor(MV), lucid_long)
m1.robust <- commarobust(m1, se_type = "CR2", clusters = 
                          na.omit(lucid_long[,c("treatment", "outcome", "mvscore","subject")])$subject)
m1.robust %>%
  modelsummary()

```

# Table F3 ----------------------------------------------------------------

```{r}
fit.2 <- lm(outcome.2 ~ treatment.2*mvscore.1 + 
             treatment.2*mvscore.2, data.frame(lucid_long) %>% 
              reshape(idvar = "subject", 
                      timevar = "round", 
                      direction = "wide"), MV.1 != 0 & MV.2 != 0)

fit.2 %>%
  modelsummary()
```

# Figure G1 ----------------------------------------------------------------

```{r}
m.fx <- meta.summaries(meta_estimates$fx, meta_estimates$se, method = "random")

data.frame(fx = c(meta_estimates$fx, m.fx$summary), 
           se = c(meta_estimates$se, m.fx$se.summary),
           study = c("Student Loans (Lucid)",
                     "KKK (Lucid)",
                     "Welfare (Lucid)",
                     "Immigration (Lucid)",
                     "KKK (Qualtrics)",
                     "Student Loans (Lucid 2)",
                     "KKK (Lucid 2)",
                     "Welfare (Lucid 2)",
                     "Immigration (Lucid 2)",
                     "Student Loans (NORC)",
                     ""),
           Sample = c(rep("Lucid", 4), "Qualtrics", 
                      rep("Lucid", 4), "NORC", ""),
           type = c(rep("Experiment", 10),
                    "Meta-Analysis Estimate")) %>%
 ggplot(aes(x = forcats::fct_reorder(study, -abs(fx)), y = fx, group = Sample)) + 
 geom_pointrange(aes(ymin = fx - 1.96*se, ymax = fx + 1.96*se)) +
 facet_grid(type ~ ., scales = "free") +
 scale_color_grey() +
 geom_hline(yintercept = 0) +
 theme_bw() + 
 coord_flip() +
 theme(text = element_text(size=15)) +
 labs(y = "DiD Estimate (Control Group SDs)", x = "") 
```

# Figure G2 ----------------------------------------------------------------

```{r}
lucid_long$mv_binary <- as.numeric(lucid_long$MV > 0)

lucid_long %>%
 mutate(MV = plyr::mapvalues(mv_binary, c(0, 1), c("No MV", "MV"))) %>%
 mutate(t = plyr::mapvalues(treatment, c(0, 1), c("Control", "Treatment"))) %>%
 mutate(exp = plyr::mapvalues(EXP, 1:4, c("Student Loans", "KKK", "Welfare", "Immigration"))) %>%
 group_by(exp, MV, t) %>%
 ggplot(aes(x = outcome, group = MV, fill = MV)) +
 geom_density(alpha = .5) +
 facet_wrap(exp ~ t, ncol = 2, nrow = 4) +
 theme(legend.position = "bottom") +
 labs(x = "Outcome", y = "Density")
```

# Figure H1 ---------------------------------------------------------------

```{r}
mv_power %>%
 mutate(study = factor(study, 
                       c("MTurk1", "Qualtrics", "NORC",
                         "MTurk2", "Lucid"))) %>%
 ggplot(aes(x = AtLeast_n_CorrectMVCs, y = t_stat)) +
 geom_text(aes(label = paste("TE=", round(beta, 2), "\n",
                             "(N=", n_size, ")", sep = "")),
           position = position_nudge(y = 3.5),
           size = 2.5) +
 geom_line() +
 geom_point() +
 facet_wrap(. ~ study) +  
 labs(x = "Subset Answering >=n MVCs Correctly") +
 theme_bw() +
 labs(y = "t-statistic") +
 scale_x_continuous(expand = expand_scale(add = .75)) +
 ylim(0, 40) +
 geom_hline(yintercept = 1.96, color = "red", linetype = 2)
```

# Table I1 ----------------------------------------------------------------

```{r}
mvc.model <- estimatr::lm_robust(outcome_stdized ~ treatment*mvc_scale, 
                                 data =  lucid_2021, 
                                 clusters = id)

imc.model <- estimatr::lm_robust(outcome_stdized ~ treatment*imc_scale, 
                                 data =  lucid_2021, 
                                 clusters = id)

modelsummary(list(mvc.model, imc.model))
```

# Figure I1 ---------------------------------------------------------------

```{r}
mvc.exps <- sapply(1:4, function(x) 
 estimatr::lm_robust(outcome_stdized ~ treatment*mvc_scale, data =  subset(lucid_2021, experiment == x),
                     clusters = id, fixed_effects = ~ round + mvc))

imc.exps <- sapply(1:4, function(x) 
 estimatr::lm_robust(outcome_stdized ~ treatment*imc_scale, data =  subset(lucid_2021, experiment == x),
                     clusters = id, fixed_effects = ~ round + mvc))

meta.coefs <- c(mvc.exps[,1]$coefficients[3] - imc.exps[,1]$coefficients[3],
                mvc.exps[,2]$coefficients[3] - imc.exps[,2]$coefficients[3],
                mvc.exps[,3]$coefficients[3] - imc.exps[,3]$coefficients[3],
                mvc.exps[,4]$coefficients[3] - imc.exps[,4]$coefficients[3])

meta.ses <- c(sqrt(mvc.exps[,1]$std.error[3]^2 + imc.exps[,1]$std.error[3]^2),
              sqrt(mvc.exps[,2]$std.error[3]^2 + imc.exps[,2]$std.error[3]^2),
              sqrt(mvc.exps[,3]$std.error[3]^2 + imc.exps[,3]$std.error[3]^2),
              sqrt(mvc.exps[,4]$std.error[3]^2 + imc.exps[,4]$std.error[3]^2))

summary(meta.summaries(meta.coefs, meta.ses, method = "random"))

tibble(effects = c(mvc.exps[,1]$coefficients[3], imc.exps[,1]$coefficients[3],
                   mvc.exps[,2]$coefficients[3], imc.exps[,2]$coefficients[3],
                   mvc.exps[,3]$coefficients[3], imc.exps[,3]$coefficients[3],
                   mvc.exps[,4]$coefficients[3], imc.exps[,4]$coefficients[3],
                   meta.summaries(meta.coefs, meta.ses, method = "random")$effects,
                   meta.summaries(meta.coefs, meta.ses, method = "random")$summary),
       se = c(mvc.exps[,1]$std.error[3], imc.exps[,1]$std.error[3],
              mvc.exps[,2]$std.error[3], imc.exps[,2]$std.error[3],
              mvc.exps[,3]$std.error[3], imc.exps[,3]$std.error[3],
              mvc.exps[,4]$std.error[3], imc.exps[,4]$std.error[3],
              meta.summaries(meta.coefs, meta.ses, method = "random")$stderrs,
              meta.summaries(meta.coefs, meta.ses, method = "random")$se.summary),
       experiment = as.factor(c(rep("Student Loans", 2),
                                rep("KKK", 2),
                                rep("Welfare", 2),
                                rep("Immigration", 2),
                                "Student Loans", "KKK", "Welfare", "Immigration", 
                                " ")),
       type = c("CATE", "CATE", "CATE", "CATE", 
                "CATE", "CATE", "CATE", "CATE", 
                "Difference in CATEs (MVC vs. IMC)", "Difference in CATEs (MVC vs. IMC)", 
                "Difference in CATEs (MVC vs. IMC)", "Difference in CATEs (MVC vs. IMC)", 
                "Random Effects Meta-Analysis"),
       Comparison = c("MVC", "IMC", "MVC", "IMC", "MVC", "IMC", "MVC", "IMC",
                      "", "", "", "", 
                      "")) %>%
 ggplot(aes(x = forcats::fct_reorder(experiment,effects), 
            y = effects, group = Comparison, shape = Comparison,
            label = Comparison)) +
 geom_hline(yintercept = 0, linetype = "dashed", alpha = .25) +
 geom_linerange(aes(ymin = effects - 1.96*se,
                    ymax = effects + 1.96*se),
                position = position_dodge(.75)) +
 geom_point(color = "black", 
            size = 5, position = position_dodge(.75)) +
 theme_ipsum_rc() +
 coord_flip() +
 facet_wrap(. ~ type, scales = "free_y", 
            nrow = 3) +
 labs(x = "", y = "Estimate") +
 scale_y_continuous(n.breaks  = 10) +
 theme(legend.position = "bottom") 
```

# Figure I2 ---------------------------------------------------------------

```{r}

# fmc comparison
fmc.mvc <- estimatr::lm_robust(fmc ~ mvc_scale, data =  lucid_2021, 
                               clusters = id, fixed_effects = ~ round + experiment + mvc)

fmc.imc <- estimatr::lm_robust(fmc ~ imc_scale, data =  lucid_2021, 
                               clusters = id, fixed_effects = ~ round + experiment + mvc)

est.diff <- (fmc.mvc$coefficients[1] - fmc.imc$coefficients[1])
se.diff <- sqrt(fmc.mvc$std.error[1]^2 + fmc.imc$std.error[1]^2)
t.stat.diff <- est.diff/se.diff

lucid_2021$log_exp_duration <- log(lucid_2021$exp_duration+.01)
lucid_2021$log_dv_duration <- log(lucid_2021$dv_duration+.01)


# exp duration comparison
exp.duration.mvc <- estimatr::lm_robust(log_exp_duration ~ mvc_scale, data =  lucid_2021, 
                                        clusters = id, fixed_effects = ~ round + experiment + mvc)

exp.duration.imc <- estimatr::lm_robust(log_exp_duration ~ imc_scale, data =  lucid_2021, 
                                        clusters = id, fixed_effects = ~ round + experiment + mvc)

est.diff <- (exp.duration.mvc$coefficients - exp.duration.imc$coefficients)
se.diff <- sqrt(exp.duration.mvc$std.error^2 + exp.duration.imc$std.error^2)
t.stat.diff <- est.diff/se.diff

# dv duration comparison
dv.duration.mvc <- estimatr::lm_robust(log_dv_duration ~ mvc_scale, data =  lucid_2021, 
                                       clusters = id, fixed_effects = ~ round + experiment + mvc)

dv.duration.imc <- estimatr::lm_robust(log_dv_duration ~ imc_scale, data =  lucid_2021, 
                                       clusters = id, fixed_effects = ~ round + experiment + mvc)

est.diff <- (dv.duration.mvc$coefficients - dv.duration.imc$coefficients)
se.diff <- sqrt(dv.duration.mvc$std.error^2 + dv.duration.imc$std.error^2)
t.stat.diff <- est.diff/se.diff

total.duration.imc <- estimatr::lm_robust(log_total_duration ~ imc_scale, lucid_2021)
total.duration.mvc <- estimatr::lm_robust(log_total_duration ~ mvc_scale, lucid_2021)

bind_rows(broom::tidy(total.duration.imc), 
          broom::tidy(total.duration.mvc), 
          broom::tidy(exp.duration.imc), 
          broom::tidy(exp.duration.mvc), 
          broom::tidy(dv.duration.imc), 
          broom::tidy(dv.duration.mvc)) %>%
 filter(term != "(Intercept)") %>%
 mutate(term = plyr::mapvalues(term, c("mvc_scale", "imc_scale"), 
                               c("MVC Score (0-1)", "IMC Score (0-1)"))) %>%
 mutate(outcome = plyr::mapvalues(outcome, c("log_total_duration",
                                             "log_exp_duration",
                                             "log_dv_duration"),
                                  c("Total Duration",
                                    "Treatment Duration", "Outcome Duration"))) %>%
 ggplot(aes(x = term, y = estimate)) +
 geom_pointrange(aes(ymin = estimate - 1.96*std.error,
                     ymax = estimate + 1.96*std.error)) +
 coord_flip() +
 facet_wrap(. ~ fct_reorder(outcome, c(1,2,3, 1,2,3)), nrow = 3) +
 labs(x = "", y = "Marginal Effect on Total Duration (Logged)") +
 theme_ipsum_rc() 

```

# Figure I3 ---------------------------------------------------------------

```{r}
bind_rows(broom::tidy(fmc.mvc), 
          broom::tidy(fmc.imc)) %>%
 mutate(term = plyr::mapvalues(term, c("mvc_scale", "imc_scale"), 
                               c("MVC Score (0-1)", "IMC Score (0-1)"))) %>%
 ggplot(aes(x = term, y = estimate)) +
 geom_pointrange(aes(ymin = estimate - 1.96*std.error,
                     ymax = estimate + 1.96*std.error)) +
 coord_flip() +
 ylim(c(0, .4)) +
 labs(x = "", y = "Marginal Effect on Factual Manipulation Check Passage") +
 theme_ipsum_rc()

```

# Figure I4 ---------------------------------------------------------------

```{r}
fit_mvc <- lm(mvc_high ~ scale(age) + scale(income) + scale(educ) + scale(ideology) + scale(polint) + scale(pid7) + scale(savviness), data = lucid2021w)

fit_imc <- lm(imc_high ~ scale(age) + scale(income) + scale(educ) + scale(ideology) + scale(polint) + scale(pid7) + scale(savviness), data = lucid2021w)

broom::tidy(fit_mvc) %>% bind_cols(outcome = "MVC High-Performers") %>%
 bind_rows(broom::tidy(fit_imc) %>% bind_cols(outcome = "IMC High-Performers")) %>%
 filter(term != "(Intercept)") %>%
 mutate(term = plyr::mapvalues(term, c("scale(savviness)", "scale(polint)", "scale(pid7)",
                                       "scale(income)", "scale(ideology)", "scale(educ)", "scale(age)"),
                               c("Savviness", "Political Interest", "Partisanship", "Income",
                                 "Ideology", "Education", "Age"))) %>%
 ggplot(aes(x = term, y = estimate, group = outcome, shape = outcome)) +
 geom_pointrange(aes(ymin = estimate - std.error*1.96,
                     ymax = estimate + std.error*1.96),
                 position = position_dodge(.5)) +
 coord_flip() +
 geom_hline(yintercept = 0, linetype = "dashed") +
 theme_ipsum_rc() +
 theme(legend.position = "bottom") +
 labs(shape = "", y = "Standardized Marginal Effect Estimate", x = "")
```

# Figure J1 ---------------------------------------------------------------

```{r echo = FALSE}
### Start 2SLS analyses
library(ivreg)
library(arm)
df <- data.frame()
latencies <- c(0, 15, 25, 35, 45)

for (latency in latencies) { 
 
 treatment <- ifelse(lucid_wide$SLexpTR1 == 1, 1, 0)
 treatment.timer <- lucid_wide$Q146_3
 treatment.complier <- ifelse(treatment == 1 & 
                               treatment.timer > latency, 1, 0)
 
 
 treatment.complier.2 <- ifelse(treatment == 1 & 
                                 lucid_wide$fmc1 == 1, 1, 0)
 
 
 control <- ifelse(lucid_wide$SLexpTR1 == 0, 1, 0)
 control.timer <- lucid_wide$Q144_3
 control.complier <- ifelse(control == 1 & 
                             control.timer > latency, 1, 0)
 control.complier.2 <- ifelse(control == 1 & 
                               lucid_wide$fmc1 == 1, 1, 0)
 
 outcome <- (lucid_wide$SLexpDV - mean(lucid_wide$SLexpDV[control == 1], na.rm = T))/sd(lucid_wide$SLexpDV[control == 1], na.rm = T)  
 
 c1 <- coef(lm(outcome ~ treatment))["treatment"]
 c2 <- coef(ivreg(outcome ~ treatment.complier | treatment))["treatment.complier"]
 c3 <- coef(lm(outcome ~ control))["control"]
 c4 <- coef(ivreg(outcome ~ control.complier | control))["control.complier"]
 c5 <- coef(ivreg(outcome ~ treatment.complier.2 | treatment))["treatment.complier.2"]
 c6 <- coef(ivreg(outcome ~ control.complier.2 | control))["control.complier.2"]
 
 s1 <- se.coef(lm(outcome ~ treatment))["treatment"]
 s2 <- sqrt(diag(vcov(ivreg(outcome ~ treatment.complier | treatment))))[2]
 s3 <- se.coef(lm(outcome ~ control))["control"]
 s4 <- sqrt(diag(vcov(ivreg(outcome ~ control.complier | control))))[2]
 s5 <- sqrt(diag(vcov(ivreg(outcome ~ control.complier.2 | control))))[2]
 s6 <- sqrt(diag(vcov(ivreg(outcome ~ treatment.complier.2 | treatment))))[2]
 
 df <- rbind(df, data.frame(timer = latency, 
                            study = "Student Loans",
                            model = c("CACE (T=1)", 
                                      "CACE (T=0)",
                                      "CACE (T=1)",
                                      "CACE (T=0)"), 
                            type = c(rep("Time", 2), rep("FMC", 2)),
                            rbind(c(c2, s2), c(c4, s4),
                                  c(c5, s5), c(c6, s6))))
 
}

df$timer <- as.factor(df$timer)

for (latency in latencies) { 
 
 treatment <- ifelse(lucid_wide$KKKexpC1 == 1, 1, 0)
 treatment.timer <- lucid_wide$Q182_3
 treatment.complier <- ifelse(treatment == 1 & 
                               treatment.timer > latency, 1, 0)
 
 treatment.complier.2 <- ifelse(treatment == 1 & 
                                 lucid_wide$fmc2 == 1, 1, 0)
 
 
 control <- ifelse(lucid_wide$KKKexpC1 == 0, 1, 0)
 control.timer <- lucid_wide$Q184_3
 control.complier <- ifelse(control == 1 & 
                             control.timer > latency, 1, 0)
 control.complier.2 <- ifelse(control == 1 & 
                               lucid_wide$fmc2 == 1, 1, 0)
 
 outcome <- (lucid_wide$KKKexpDV - mean(lucid_wide$KKKexpDV[control == 1], na.rm = T))/sd(lucid_wide$KKKexpDV[control == 1], na.rm = T)
 
 c1 <- coef(lm(outcome ~ treatment))["treatment"]
 c2 <- coef(ivreg(outcome ~ treatment.complier | treatment))["treatment.complier"]
 c3 <- coef(lm(outcome ~ control))["control"]
 c4 <- coef(ivreg(outcome ~ control.complier | control))["control.complier"]
 c5 <- coef(ivreg(outcome ~ treatment.complier.2 | treatment))["treatment.complier.2"]
 c6 <- coef(ivreg(outcome ~ control.complier.2 | control))["control.complier.2"]
 
 s1 <- se.coef(lm(outcome ~ treatment))["treatment"]
 s2 <- sqrt(diag(vcov(ivreg(outcome ~ treatment.complier | treatment))))[2]
 s3 <- se.coef(lm(outcome ~ control))["control"]
 s4 <- sqrt(diag(vcov(ivreg(outcome ~ control.complier | control))))[2]
 s5 <- sqrt(diag(vcov(ivreg(outcome ~ control.complier.2 | control))))[2]
 s6 <- sqrt(diag(vcov(ivreg(outcome ~ treatment.complier.2 | treatment))))[2]
 
 df <- rbind(df, data.frame(timer = latency, 
                            study = "KKK",
                            model = c("CACE (T=1)", 
                                      "CACE (T=0)",
                                      "CACE (T=1)",
                                      "CACE (T=0)"), 
                            type = c(rep("Time", 2), rep("FMC", 2)),
                            rbind(c(c2, s2), c(c4, s4),
                                  c(c5, s5), c(c6, s6))))
 
}

df$timer <- as.factor(df$timer)

###

for (latency in latencies) { 
 
 treatment <- ifelse(lucid_wide$APexpLazy == 0, 1, 0)
 treatment.timer <- lucid_wide$Q54_3
 treatment.complier <- ifelse(treatment == 1 & 
                               treatment.timer > latency, 1, 0)
 
 treatment.complier.2 <- ifelse(treatment == 1 & 
                                 lucid_wide$fmc3 == 1, 1, 0)
 
 
 control <- ifelse(lucid_wide$APexpLazy == 1, 1, 0)
 control.timer <- lucid_wide$Q12_3
 control.complier <- ifelse(control == 1 & 
                             control.timer > latency, 1, 0)
 control.complier.2 <- ifelse(control == 1 & 
                               lucid_wide$fmc3 == 1, 1, 0)
 
 outcome <- (lucid_wide$APexpDV - mean(lucid_wide$APexpDV[control == 1], na.rm = T))/sd(lucid_wide$APexpDV[control == 1], na.rm = T)
 
 c1 <- coef(lm(outcome ~ treatment))["treatment"]
 c2 <- coef(ivreg(outcome ~ treatment.complier | treatment))["treatment.complier"]
 c3 <- coef(lm(outcome ~ control))["control"]
 c4 <- coef(ivreg(outcome ~ control.complier | control))["control.complier"]
 c5 <- coef(ivreg(outcome ~ treatment.complier.2 | treatment))["treatment.complier.2"]
 c6 <- coef(ivreg(outcome ~ control.complier.2 | control))["control.complier.2"]
 
 s1 <- se.coef(lm(outcome ~ treatment))["treatment"]
 s2 <- sqrt(diag(vcov(ivreg(outcome ~ treatment.complier | treatment))))[2]
 s3 <- se.coef(lm(outcome ~ control))["control"]
 s4 <- sqrt(diag(vcov(ivreg(outcome ~ control.complier | control))))[2]
 s5 <- sqrt(diag(vcov(ivreg(outcome ~ control.complier.2 | control))))[2]
 s6 <- sqrt(diag(vcov(ivreg(outcome ~ treatment.complier.2 | treatment))))[2]
 
 df <- rbind(df, data.frame(timer = latency, 
                            study = "Welfare",
                            model = c("CACE (T=1)", 
                                      "CACE (T=0)",
                                      "CACE (T=1)",
                                      "CACE (T=0)"), 
                            type = c(rep("Time", 2), rep("FMC", 2)),
                            rbind(c(c2, s2), c(c4, s4),
                                  c(c5, s5), c(c6, s6))))
 
}

for (latency in latencies) { 
 
 treatment <- ifelse(lucid_wide$ImmExpC2 == 0, 1, 0)
 treatment.timer <- lucid_wide$Q160_3
 treatment.complier <- ifelse(treatment == 1 & 
                               treatment.timer > latency, 1, 0)
 
 treatment.complier.2 <- ifelse(treatment == 1 & 
                                 lucid_wide$fmc4 == 1, 1, 0)
 
 
 control <- ifelse(lucid_wide$ImmExpC2 == 1, 1, 0)
 control.timer <- lucid_wide$Q162_3
 control.complier <- ifelse(control == 1 & 
                             control.timer > latency, 1, 0)
 control.complier.2 <- ifelse(control == 1 & 
                               lucid_wide$fmc4 == 1, 1, 0)
 
 outcome <- (lucid_wide$ImmEXPDVscale - mean(lucid_wide$ImmEXPDVscale[control == 1], na.rm = T))/sd(lucid_wide$ImmEXPDVscale[control == 1], na.rm = T)
 
 c1 <- coef(lm(outcome ~ treatment))["treatment"]
 c2 <- coef(ivreg(outcome ~ treatment.complier | treatment))["treatment.complier"]
 c3 <- coef(lm(outcome ~ control))["control"]
 c4 <- coef(ivreg(outcome ~ control.complier | control))["control.complier"]
 c5 <- coef(ivreg(outcome ~ treatment.complier.2 | treatment))["treatment.complier.2"]
 c6 <- coef(ivreg(outcome ~ control.complier.2 | control))["control.complier.2"]
 
 s1 <- se.coef(lm(outcome ~ treatment))["treatment"]
 s2 <- sqrt(diag(vcov(ivreg(outcome ~ treatment.complier | treatment))))[2]
 s4 <- sqrt(diag(vcov(ivreg(outcome ~ control.complier | control))))[2]
 s5 <- sqrt(diag(vcov(ivreg(outcome ~ control.complier.2 | control))))[2]
 s6 <- sqrt(diag(vcov(ivreg(outcome ~ treatment.complier.2 | treatment))))[2]
 
 df <- rbind(df, data.frame(timer = latency, 
                            study = "Immigration",
                            model = c("CACE (T=1)", 
                                      "CACE (T=0)",
                                      "CACE (T=1)",
                                      "CACE (T=0)"), 
                            type = c(rep("Time", 2), rep("FMC", 2)),
                            rbind(c(c2, s2), c(c4, s4),
                                  c(c5, s5), c(c6, s6))))
 
}

df$condition <- c(rep(c("CACE (Treatment)", 
                        "CACE (Control)"),
                      10),
                  rep(c("CACE (Free Speech)", 
                        "CACE (Safety)"),
                      10),
                  rep(c("CACE (Unlucky)", 
                        "CACE (Lazy)"),
                      10),
                  rep(c("CACE (High skilled)", 
                        "CACE (Low skilled)"),
                      10))

df$timer <- as.factor(df$timer)
df$timer <- plyr::mapvalues(df$timer, levels(df$timer), c("Full Sample", "15", "25", "35", "45"))

df %>% filter(type == "FMC" & timer == "Full Sample" | type == "Time") %>% 
 mutate(treatment = abs(treatment.complier),
        se = treatment.complier.1) %>%
 ggplot(aes(x = condition, y = treatment, color = timer, shape = type)) +
 geom_pointrange(aes(ymin = treatment - 1.96*se,
                     ymax = treatment + 1.96*se), size = .5, 
                 position = position_dodge(.5)) +
 theme_bw() +
 theme(axis.text.x = element_text(angle = 90)) +
 facet_wrap(. ~ study, scales = "free") +
 scale_color_grey(start = .2, end = .8) +
 labs(x = "Complier Subgroup", y = "Absolute Effect Estimate (SDs)",
      color = "Latency Threshold",
      shape = "Attentiveness Measure") +
 theme(legend.position = "bottom") +
 guides(colour = guide_legend(override.aes = list(shape = 17)))

```

# Figure K1 ---------------------------------------------------------------

```{r echo = FALSE}
library(interflex)

interflex(estimator = "binning",
          Y = "outcome",
          D = "treatment",
          X = "mvscore",
          data = data.frame(lucid_long),
          na.rm = T,
          nbins = 2)

interflex(estimator = "kernel",
          Y = "outcome",
          D = "treatment",
          X = "mvscore",
          bw = .40,
          data = data.frame(lucid_long),
          na.rm = T)
```
